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A new technique is presented for modifying the Einstein evolution equations off the constraint 
hypersurface. With this approach the evolution equations for the constraints can be specified freely. 
The equations of motion for the gravitational field variables are modified by the addition of terms 
that are linear and nonlocal in the constraints. These terms are obtained from solutions of the 
linearized Einstein constraints. 



The Einstein equations separate into a set of evolu- 
tion equations and a set of constraints. The evolution 
equations are partial differential equations (PDE's) that 
determine how the gravitational field variables g a b (the 
spatial metric) and K a b (the extrinsic curvature) evolve 
forward in time. The constraint equations are PDE's 
that the field variables must satisfy at each instant of 
time. From a Hamiltonian point of view, the evolution 
equations define solution trajectories in phase space with 
coordinates g a b and momenta K a b. Physical trajectories 
are those that lie in the constraint hypersurface, or sub- 
space, of the gravitational phase space. 

Einstein's theory of gravity is a "first class" theory, 
that is, the time derivatives of the constraints are linear 
combinations of the constraints. This property implies 
that, analytically, the constraints will hold at each in- 
stant of time if they hold at the initial time. However, for 
numerically generated solutions of the theory the initial 
data will not satisfy the constraints precisely and numer- 
ical errors will kick the phase space trajectory away from 
the constraint hypersurface. This is a critical problem 
for numerical modeling because the Einstein evolution 
equations, as they are usually written, admit solutions 
that rapidly diverge away from the constraint hypersur- 
face Any numerical scheme that evolves the grav- 
itational field data using the evolution equations in one 
of their traditional forms will eventually fail to produce 
physically meaningful results. Inevitably the numerical 
solution will choose to follow a trajectory that violates 
the constraints. 

A number of strategies have been devised to ad- 
dress this problem. One approach is to modify the 
theory off the constraint hypersurface by adding linear 
combinations of constraints to the evolution equations 
|H IE IE IE El In this way one hopes to alter the so- 
lution trajectories so that they are better behaved away 
from the constraint hypersurface. We will use the termi- 
nology "off-shell" to refer to solution trajectories that lie 
off the constraint hypersurface. 

The strategy discussed in this paper is of this sort. 
We add terms proportional to the constraints to the Ein- 
stein evolution equations in such a way that the evolution 
equations for the constraints can be freely specified. In 
principle we can eliminate all constraint violating modes 
by demanding, for example, that the time derivatives of 



the constraints should vanish. The price we pay for this 
degree of control over the unphysical, off shell solutions 
is that the terms added to the evolution equations are 
nonlocal. They are determined through the solution of 
an elliptic system of PDE's. 

Another strategy for keeping a numerically generated 
solution from diverging away from the constraint hyper- 
surface is constrained evolution. In this scheme the con- 
straints are used in place of certain evolution equations to 
update some of the gravitational field variables in time. 
This approach has worked well for spherically and ax- 
isymmetric problems [E IE ll E llEl- A closely related 
idea is constraint projection yj ll3|, Ll| • With constraint 
projection one evolves the full set of field variables using 
the evolution equations, then periodically (perhaps every 
timestep) solves the constraints to project the solution 
back to the constraint hypersurface. Both constrained 
evolution and constraint projection require the solution 
of the constraint equations during the course of evolution. 
For these approaches to be viable, the constraints must 
be expressed as an elliptic system of PDE's. From a com- 
putational perspective, our strategy is closely related to 
constraint projection since we also solve an elliptic sys- 
tem of PDE's at every (or nearly every) timestep. In 
fact, the PDE's that we solve are the linearized Einstein 
constraints. 

It will be useful to give an overview of our procedure for 
modifying the off-shell solutions in a formal, general con- 
text. Consider a theory, like general relativity, described 
by a set of first class constraints Ca- Let tp^ denote the 
basic field variables. These variables satisfy first order 
in time differential equations of motion, ip^ = (?/v)oid rhs> 
where the "old right-hand sides" (^)oidrhs are functions 
of and their spatial derivatives. We have included the 
descriptor "old" since we will soon create "new" right- 
hand sides by adding functions of the constraints. The 
evolution equations for the constraints are obtained from 
the evolution equations for the ^'s by differentiating the 
constraints in time. This yields Ca = (Cyt)oidrhs where 
the right-hand sides (rhs's) are given by 

(d-Ooldrhs = T-^-(Vv)old rhs ■ (1) 

The expression SCa/SiPh is the Frechet derivative of the 
constraints with respect to the field variables 01 • It 
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satisfies Ca{$ + f) — Ca(VO = (SCa/ 'fo/v)°V" in the limit 
as the norm of c M goes to zero, (ip^ and cr^ are defined 
as vectors in a suitable Banach space.) If Ca depends on 
spatial derivatives of ip^, then the Frechet derivative is a 
differential operator. 

Now let us split the basic variables ip^ into two sets, 
<pA and Xi- Note that there are as many </>'s as there are 
constraints. With this splitting, Eq. Q becomes 

(d-Ooldrhs = -rT^"(0s)oldrhs + (xO°ld rhs ■ (2) 

<5</>B OXi 

Next, we replace the old equations of motion for the (f>'s 
with new equations of motion. This leads to new equa- 
tions of motion for the constraints, 

(OOncwrhs = FT (4>B )ncw rhs + "T (Xi)oldrhs • (3) 

0(pB OXi 

By subtracting the previous two results, we find 



AA = TT— *B 

0<pB 



(4) 



where $,4 is the difference between the new and old rhs's 
for the variables (p A , ®a = (^u)newrhs - (0A)oidrhs, and 
Aa is the difference between the new and old rhs's for 
the constraints, = (CA)newrhs - (Cyt)oidrhs- In terms 
of the original field variables ip^, the new equations of 
motion are 
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(5) 



Now we turn the reasoning around. We do not actu- 
ally choose new equations of motion for the 4>'s. Instead, 
we specify new evolution equations for the constraints 
by freely choosing the expressions (Ca) new rhs- The func- 
tions Aa are then determined, and Eqs. are solved for 
The new equations of motion for the original field 
variables are given by Eqs. JSJ. 

Because the theory is first class, (CA)oid rhs is a linear 
combination of constraints. Let us choose (C J 4) ne wrhs to 
be a linear combination of constraints as well. Then Aa 
is a linear combination of constraints and, according to 
Eq. (0} , &a is a (possibly nonlocal) linear combination of 
constraints. It follows that the new equations of motion 
for ipn differ from the old equations by a linear combina- 
tion of constraints. 

Equations I0J are the linearized constraints. To be 
precise, consider a field configuration ip^ that does not 
satisfy the constraints and let -tp^ = ip^ — ^ n- If Vy satis- 
fies the constraints then to linear order, ^ satisfies Ca — 
(SCa/SiP^^. This is Eq. © with tf„ = (S^/S<p A )^A 
and Ca = A a- 

With the procedure outlined above we can freely spec- 
ify the rhs's of the constraint evolution equations, as long 
as they are linear combination of constraints. We then 
solve Eq. (0J for <&a and modify the equations of motion 



for the t/>'s to Eq. JSJ. In this way we leave the equations 
of motion for the basic field variables unchanged on the 
constraint hypersurface, but we can modify their off-shell 
form to eliminate the constraint violating modes. 

Let us apply this formalism to general relativity. The 
basic field variables (the ip's) are the spatial metric g a b 
and the extrinsic curvature K a b- The equations of motion 
as written by York [l^l are d±g a b = (d±g a b)o\d rhs and 
d±K a .b = (dxK a b)aidrha, where 

(dLSah) old rhs = -2aK ab , (6a) 
{d±Kab) oldlhs = a(KK ab ~2K ac K^ + R ab ) 

-D a D b a. (6b) 

Here, a is the lapse function, D a is the spatial covariant 
derivative, and R a b is the spatial Ricci tensor. The op- 
erator d± = dt — Cp is the difference between the time 
derivative and the Lie derivative along the shift vector 
/3 a . This operator plays the role of the "dot" in the for- 
mal analysis. 

The constraints for general relativity are 



H = K 2 -K ab K ab + R 



M a = D b K b a - D a K . 
With the equations of motion we find 



(7a) 
(7b) 



(<9iW)oldrh s = 2aKH-2aD a M a -4M a D a a , (8a) 
{d ± M a )oU rhs = aKM a - HD a a - aD a H/2 . (8b) 

for the rhs's of the constraint evolution equations. 

In order to proceed, we must select a subset of the vari- 
ables g a b, K a b to play the role of the cp's. We will use the 
conformal transverse-traceless decomposition developed 
by Lichnerowicz and York for solving the initial value 
problem (see, for example, Ref. 0])- To begin, we split 
the metric and extrinsic curvature into 



9ab 



<p 4 g a b 



K ab = <P A ab + -if g ah r 



(9b) 



where tp is the conformal factor, g a b is the conformal 
metric, r is the trace of the extrinsic curvature, and A a b 
is symmetric and trace free. Note that these definitions 
are invariant under the conformal transformation [Isl fl9| 
g a b -> £, 4 g a b, y —> £~V> A ab -> £,~ 2 A ab , r -> t. 

The tensor A a b can be decomposed in terms of a sym- 
metric, transverse, traceless tensor B a b and a vector w a , 



Aab = {lLw) ab + B ab 



(10) 



The operator IL is defined by (Kw) ab = D a w b + DbW a — 
2g a bD c w c /3. It follows that w a satisfies the elliptic equa- 
tion D b (ILw)ab = D h A a b- The conformal transformation 
rule for B ab is B ab — ► i~ 2 B ab and the rule for w a is de- 
fined by the relation (]Lw) a b — > £, 2 (]Lw)ab- 
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Let tp and w a play the role of the (fi's in Eqs. J2J and JSJ. 
The Frechet derivatives that appear in those equations 
are somewhat tedious but straightforward to compute. 
If we let $ and W a denote the unknowns (the <&a's) in 
Eq. we find the following results: 

A = -8D 2 ($/tp)-2K ab (EW) ab /tp 2 



UK 2 



Aa = D b 



12K ab K ab 

2 



4R] $/V , (11a) 



(EW) ab /p 



-6 (D b K ab - D a K/3) $/tp . 



(lib) 



Here, Ao and A a are the differences between the new and 
old rhs's of the evolution equations for the constraints: 

Ao = (d ± H) newihs -(d ± H) oldlhs , (12a) 

A a EE (d±M a ) newlhs - (d±M a ) o}dlhs ■ (12b) 

The unknowns in Eqs. (|llfl are the differences between 
the new and old rhs's for the tp and w a equations of mo- 
tion. From the formal Eq. JSJ , we find the new equations 
of motion 



(dj_g a b) now r hs 

{d±K ab ) nev/ rhs 



(d±g ab ) old rhs + ^Jab^/p , (13a) 
(d^K ab ) oldlhs + (EW) ab /p 2 
-2{K ab - Kg ab )<5>/p (13b) 



for the metric and extrinsic curvature. 

Let us assume that the new equations of motion for tp 
and w a , like the old, are conformally invariant. Then we 
see that $ and W a inherit the conformal transformation 
properties $ -> and {EW) ab -> £,- 2 {EW) ab . It 

follows that Eqs. and (|13|l are invariant under con- 
formal transformations. In other words, these equations 
hold for any choice of splitting of the physical metric 
g ab into a conformal factor tp A and conformal metric g ab . 
For simplicity we can choose the conformal factor to be 
unity, tp = 1, and the conformal metric to coincide with 
the physical metric, g ab = g ab . Then Eqs. (|11|) become 

A = -8L> 2 $ - 2K ab (EW) ab 

- [AK 2 - \2K ab K ab + 4R] $ , (14a) 
A a = D b (EW) ab - 6 (D b K ab - D a K/3) $ . (14b) 

This is an elliptic system of linear PDE's for $ and W a . 
The new equations of motion l|13l) become 



(d±g a b, 
(d±K abi 



new rhs 
new rhs 



{d±9ab) oldlhs + 4.g Q& $ , (15a) 
(9±K ab ) oldths + (EW) ab 
-2(K ab - Kg ab )$ . (15b) 



This is our main result. We can now specify new rhs's for 
the Hamiltonian and momentum constraint equations of 
motion. The A's are found from Eqs. (HJ, and Eqs. 1(14^1 
are solved for $ and W a . These results are used in the 
new equations of motion, Eqs. i|15|) . Note that these 



equations contain only the physical metric and extrin- 
sic curvature — all references to the conformal transverse- 
traceless splitting have disappeared. 

What follows is a numerical demonstration that 
Eqs. i|14fl. tEH allow us the freedom to prescribe the 
evolution of the constraints by altering the Einstein equa- 
tions off-shell. For this test we use periodic identification 
in the x, y, and z coordinate directions with periods of 
2ir. Thus, for now, we intentionally avoid facing the very 
important issue of boundary conditions. We use initial 
data that violates the constraints. The spatial metric is 
flat with Cartesian coordinates, g ab — S ab . The diago- 
nal components of the extrinsic curvature are given by 
K xx = K yy = K zz = A/3. The off-diagonal component 



K 



xy 



Si cos 2 (z) + e 2 cos(a;) cos(y) 



(16) 



The components K yz and K zx are obtained from K xy by 
cyclic permutations of x, y and z. The initial data is 
evolved with unit lapse a — 1 and vanishing shift (3 a — 
0. For the test results shown here we use the values 
A = 0.02, £i = 0.01, and e 2 = 0.0005. 

Figure^ shows the common logarithm of the L2 norm 
of the constraints as a function of time, with and without 
the off-shell modification terms in Eqs. Ijl5|l . For the new 
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FIG. 1: Common logarithm of the L2 norms of the Hamil- 
tonian and momentum constraints, using the old and new 
equations of motion. 



constraint equations we have chosen 

-0AH - L P H 



(d±M a 



new rhs 
new rhs 



(17a) 

0AM a - CpM a , (17b) 



so that at each point in space H. and M a have time de- 
pendence ~ exp(— 0.4t). The Li norms are defined with- 
out any factors of g ab or yfg. That is, we define L2H = 
VE^ 2 /^ 3 where the sum extends over the iV 3 grid 



4 



points. Similarly, we define L 2 M = M a M a /N 3 
and include a sum over the index a. With these defini- 
tions, the only time dependence that should appear in the 
L 2 norms is the exponential decay ~ exp(— OAt). This is 
precisely what we see in Fig. with the new equations 
of motion. For the old equations of motion, we see strong 
oscillations on a short time scale and exponential growth 
on a longer time scale. 

Our numerical code uses pseudospectral collocation 
[joj ] with a Fourier basis in each of the coordinate di- 
rections. Fourth — order Runge-Kutta is used for time 
stepping. The elliptic equations <|14(l are solved with the 
iterative method GMRES [2]]]. We use a left precondi- 
tioner consisting of the inverse of the diagonal part of 
the elliptic operator. One of the numerical issues that 
we face is spectral blocking [2(| . This is the phenomenon 
in which aliasing causes an unphysical increase in power 
in the highest wave number modes that can be supported 
on the grid. Filtering can help alleviate this problem. For 
the simulations shown in Fig. (JJJ , we use N = 20 colloca- 
tion points in each dimension and a filter that sets the two 
highest frequencies to zero at the end of each timestep. 
The timestep is 0.04, compared to a light-crossing time 
of approximately 2n. 

The results displayed in Fig. show that we have 
indeed modified the equations of motion off-shell in such 
a way that unwanted growth in the constraints is elimi- 
nated. Ultimately, what we would like to show is the abil- 
ity to prevent constraint growth in the first place. Our 
preliminary attempts to demonstrate this ability have not 
been completely successful, for reasons that we suspect 
are purely numerical. Although we cannot rule out the 
possibility that the combined system Eqs. i|14fl . i|15[l is 
mathematically ill-defined in some sense, the problems 
that we have encountered appear to be caused by nu- 
merical issues. One issue is spectral blocking, mentioned 
above. Another issue is the failure of our elliptic solver 
to converge to a solution under circumstances that we 
do not yet understand. We suspect that a better pre- 
conditioner will make our elliptic solver more robust and 
dependable. 

For the simulation shown in Fig. ifljl. with the new 
equations of motion, the constraints continue to drop ex- 
ponentially untile « 15. Beyond this time the constraints 
begin to suffer from high wave number variations whose 
growth counteracts the exponential drop of the longer 
wavelength modes. This breakdown is sensitive to the 
spatial resolution and amount of filtering and appears 
to be related to spectral blocking. With the old equa- 



tions of motion, the Hamiltonian and momentum con- 
straints continue to grow exponentially until t « 45. At 
that time L2H has a value of ~ lO^ 1 and the simula- 
tion breaks down. Again, this appears to be related to 
spectral blocking. We plan to study these issues in more 
detail and present further numerical tests in a future pub- 
lication. 
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